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Abstract 

In this paper, we developed a new parametrization method to calculate the localization length in one 
dimensional Anderson model with diagonal disorder. This method can avoid the divergence difficulty en 
countered in the conventional methods and significantly save computing time as well. 



* colinkk@pku.edu.cn 



1 



L INTRODUCTION 

The tight-binding Anderson model HI] is a standard model for disordered systems. It predicted 
that with random on-site energies whose strength of disorder is above a critical value, the electrons 
are localized in a certain spatial region. The transition from metal to insulator by increasing 
the strength of disorder is called Anderson transition. According to the single-parameter scaling 
theory B2D, there is no metallic regime in the disordered system with dimension d < 2 yi |4fl. 
With long-range correlated disorder, even in one-dimensional systems, there can exist extended 
electronic states and a true Anderson transition with mobility edges was observed [5|]. 

One-dimensional systems are the simplest case to study. There were a number of review articles 

Recently, articles about the 
120, ensemble-averaged 



on the problem of localization in one-dimensional systems 
general analytical methods J 1011. anomalous behaviors at band center flllL 
conductance fluctuations [130, and discrete Anderson nonlinear Schrodinger equations H14H were 
appeared in the major journals, indicating the topic is still a hot one. 

One of the important one and quasi-one dimensional real systems is the DNA molecules. The 
topics of charge transport in DNA and the feasibility of constructing DNA-based devices, have 
kindled a heated debate within the scientific community Ill5ll . Various theoretical models about 
charge transport in DNA were generalized from the original one-dimensional tight-binding model 
of Anderson yj] 



(1) 



where e, is the zth on-site energy given by a random distribution, and t# is the hopping energy from 
the ith site to the jth site. The one-dimensional Anderson model with diagonal disorder (all tp 
equal to t) expressed in terms of the Schrodinger equation is 



iffi-i + ifr M = (E - edif/i, 



(2) 



where ^, is the wave function on the ith site, E is the eigenenergy. The hopping energy t has been 
set as energy unit. Generalizations to more realistic models of DNA molecules take into account 
of correlations, random hopping energies, coupled multichains, and so on 



Of the many approaches which have been developed for numerical simu 



systems, the transfer matrix method has proved the most productive 1124 12511. Introducing a two 



ations of disordered 



component vector = (</r ; t/r,_i) f , Eq. © can be written as 

= 



E-6i -1 

1 



*Ai 1 , 



= T Y 



(3) 



where T, is the so-called transfer matrix. Here only the diagonal disorder is considered, in 
which case the transfer matrix is symplectic. With transfer matrix we can get the wavefunc- 
tion if/ L on site-L propagating from site- 1 with wavefunction tfr\ for arbitrary L in principle by 
¥ L = T i T L _ 1 ---T 1 Y 1 . 

Traditional transfer matrix method has to calculate detailed physical properties for a long chain, 
then to use scaling technique to reveal the properties under the thermodynamic limit; and it requires 
a stabilization procedure to overcome the overflow problem, typically about every twelve iteration 



j24l] . hi this paper we propose a parametrization method to deal with the transfer matrix of one- 
dimensional Anderson model with uncorrelated diagonal disorder. With this method, we directly 
calculate in the localization regime under the thermodynamic limit; and it is easy to calculate the 
localization length for arbitrary strength of disorder. Particularly, it is easier to calculate moderate 
disorder than both weak and strong disorder, in which cases it needs more numerical techniques 
to obtain the results in our method. 



II. PARAMETRIZATION OF THE TRANSFER MATRIX 



The transfer matrix T, in Eq. (O is a symplectic matrix, which satisfies the condition A'CIA = 
fl, where A' is the transpose of a In x In matrix A, and £1 is a fixed nonsingular, skew-symmetric 
matrix. Typically il is chosen to be a block matrix 



-i 



(4) 



where (D and 1 are n x n zero matrix and unit matrix respectively. It is easily shown from the 
definition that the transpose of a symplectic matrix is also a symplectic matrix. 

Let M L = T!T 2 • • • T L and v,- = E - 6 ( . Since the symplectic matrices form a group, the product 
of and M £ is symplectic and real symmetric. Thus, it can be diagonalized by an orthogonal 
matrix U(# £ ), 



U(# L )M L M L U(-# L ) = 



(5) 



where 



cos # L - sin 6 L 
k sin 9 L cos 6 L 



(6) 



The 2x2 symplectic matrix has a property that the two eigenvalues are reciprocal to each other. 
Hence, M^M L can be expressed in terms of A L and 6 L , 

cos 26 L - sin 26 L 

M^M/_ = cosh AiX + sinh Al , (7) 

^ - sin 26 L - cos 26 L 

here 1 is a 2 x 2 unit matrix and the two parameters Al and 6l play important roles to parameterize 
the transfer matrices. 

Making use of the definition of Ml, 

M^ +1 M L+1 = 

{ -i 

and substituting the expression for M^M £ in Eq. ©, we can easily obtain the recursion relations 
for A and 6, 



M'M L 



Vl+i -1 
1 



(8) 



cosh A L+ i = (1 + -^-) cosh i L + (-^- cos 2# L - v L+ \ sin 2# L ) sinh /} L 



L+l 



sinh/l L+ i cos 20/.+ 



2 2 

[(-y 1 - 1) cos 26 L - v L+ i sin 26> L ] sinh A L + cosh A L 
sinh/l L+ i sin2# L+ i = v^+i cosh/t L + (vl+i cos 20/, - sin2# L ) sinh/l L . 



(9) 

(10) 
(11) 



In the localization regime, the localization length is finite. For sufficiently long chains (de- 
pending on the strength of localization), the exponent A of the eigenvalue will be approaching to 
infinity. Hence after dividing Eq. (fTTI) by Eq. (flOl) and taking the limit A — » oo or tanh A L — > 1, we 
get the recursion relation of 20 in the localization regime, 

vl+i(1 + cos 26 l) - sin2# L 



tan 20, 



L+l 



v 2 T , (1 + cos 2Q L )/2 - v L +\ sin 26 L - cos 2# L 



(12) 



Using the basic relations of trigonometric functions, it is easy to derive a much simpler recursion 
relation of 9, Eq. CCS), from Eq. (fT2~l) . 

tan# L+1 = 1 —— . (13) 

v L+ i - tan# L 

From Eq. (fT3l) it is apparent that if vl+i = 0, = 6l± n/2, which is also indicated in Eq. (PT21) . 
Eq. (fT3l can also be written as a continued fraction 

1 



In this form we can see that, for each specified chain, 6 L is completely determined by the sequence 
of random on-site energies and the eigenenergy E in the Schrodinger equation. We should bear 
in mind that all the results we have in this paper are in the localization regime under the ther- 
modynamic limit. Therefore, although we label the first site as site 1, it does not mean that site 1 
starts from the very beginning, there are sufficiently many sites before it. Furthermore, as an initial 
input, the effect of 6\ will vanish for sufficiently large L. 

Now we derive the probability density of 6 from the known random distribution of e (or v = 
E - e). And we can calculate the localization length directly through the probability density of 9. 
Define t = tan 8. From Eq. (TT3T ), we have an integral equation for the probability density function 
of t 



J /"»CO /-*oo J 

Pc(-)= p t (t')p v (v)6[--(v-t')Wdv 

' \J —CO %J —CO * 



f 

*J — o 



Pt (t')p v (- t + t'W, (15) 



where p t (f), Pc(—\ an d Pv(v) are the probability density functions for tan#, cot#, and v respec- 
tively. With the relations of the probability density functions, 

P(&) = -\- Q Pt(i) = -A-Pcib, (16) 
cos 2 6 sin 2 6 t 

the integral equation of the probability density function p{6) becomes 

p{0) = f ' p{ff)p v {-^— + tan W. (17) 

sin 2 6 J-n/2 tan 6 

If we can obtain the solution of p(9) from Eq. ([T71) . the inverse localization length is given by 

Eq. © 

^ = r = \( a l+i ~ a l) 

'Jt/2 fOo 

p(8)p v (v)ln(\ + v 2 cos 2 9- vsin20)dvde, (18) 

12 J-oo 



L J-n/2 J-t 



where y is the so-called Lyapunov exponent. The second line is obtained when A — » oo, hence 
cosh A L+l / cosh A L —* exp(A L+l - A L ) and tanh/l L — » 1. For uncorrelated disorder considered in 
this paper, p(6) and p v (v) are independent of each other, thus the above equation is correct for 
sufficiently long chains in the localization regime. In [10], a similar expression was obtained from 
a different starting point for p v (v) as a discrete distribution and the parameter 6 was chosen to be 
real at the beginning. 



It should be mentioned that from the original Schrodinger equation we have 

zl+i = , (19) 

Vl+i - z L 

where zl = 'Al/'Al+i- This is exactly the form of Eq. ( IT3T ). The difference is that z is a complex 
number while tan 9 is real. We found that for sufficiently large L with the same random sequence, 
the real part of zl is the same as tan 9 L and its imaginary part goes to zero in the localization regime. 
Zl and tan# L are essentially the same quantity when L is sufficiently large in the localization 
regime. Thus the phase difference between ij/ L and <A i+1 is or n and the ratio of their amplitudes 
is Re((A L /t/r L+1 ) — > tan Q L for sufficiently large L in the localization regime. And the expression for 
the localization length of z in flg] can be rewritten in our symbols as the following, 



y{E) = -\ d#p(#)ln|tan<9|, (20) 

once we solve the probability density p{6) from Eq. (fTTT) . where the dependence on the eigenenergy 
E is implicitly included in p(9). The z defined in flg] is the inverse of ours, thus we have a minus 
sign in our equation. Eq. (fT~8l) and Eq. (1201) are equivalent 



p(6)p v (v) ln(l + v 2 cos 2 6 - v sin 2#)dvd# 



= \ij P( e ) ln cos2 m + JJJ dedvd cot d'p(d)Pv(v) ln[l + (v - tan 6) 2 ]S[cot 6' - (v - tan 6)] 
= -J p(#)ln|tan6>|d<9. 

However, Eq. (TT8T) is more appropriate for numerical calculation because it requires less accurate 
p{6) than Eq. (f20]). 

From the above description, we can experience benefits of our method. We do not have to 
use the recursive relation Eq. (fT3l) to calculate the sequence 6 L , nor A L . Although the recursive 
relation is very simple, the time consumption of the computation of 6 L is still the same as in 
the conventional transfer matrix method within the same accuracy. Moreover, the conventional 
transfer matrix method has a problem that it has to care about the overflow problem especially 
for strong disorder. This is not a problem in our method. We can calculate arbitrary strength 
disorder. More importantly, we can easily calculate moderate strength disorder, which is difficult 
for analytical methods. 



III. NUMERICAL TECHNIQUE 



Eq. (fT71) shall be solved numerically in a discrete matrix form by using self-consistent iterative 
procedure. In principle, we can get the probability density p(6) for any known random distribution 
p v (v). In this paper we only considered two distributions: Lorentzian and Gaussian. For the 
Lorentzian distribution, there are exact analytical results Jg] and for the Gaussian distribution, 
there are analytical results in the weak and strong disorder limits [26]. We will compare our 
numerical results with theirs. 

From the r.h.s of Eq. (fTTT ), we can see that there is a removable singularity at 6 = in both 
of the Lorentzian and Gaussian distributions in the numerical computation. Thus we need to 
use interpolation to obtain p(9) near 8 = 0. For the Lorentzian distribution we use six-points 
interpolation and for the Gaussian distribution we use seven-points interpolation with a condition 
p(0) = [pin/2) + p(-n/2)]/2. This condition can be easily derived from Eq. (fl5l>. Let t" = f + -, 
then Eq. (fT5l) becomes, 

Pc(b= f p t (t"--) Pv (t")dt". (21) 

' kJ —CO & 

For in «: \l/t\, p t {t" - l/t) « p,(-l/t), so Eq. (I2TT) can become approximately, 

Pc(-)*Pt(~) f PvV'W, (22) 

' ' *J —CO 

and we have lim p r (—) = lim pA — ). This relation can be rewritten in terms of the probabil- 
ity density of 6 as p{6 —* 0*) = p(+n/2) by using Eq. (fT6l) . We should emphasize that for the 
Lorentzian distribution this relation is not valid. The range of 6 need to use interpolation is de- 
termined by |0| < yJmn/Na, where m is a number we choose to ensure that the profile of p{6) 
near 6 = is smooth (we assume this is true for nonsingular distributions p v (v)), N is the number 

of equally spaced abscissas points where the integrands are evaluated and <x is the parameter of 

1 ?? 1 O" 

the distributions, p v (v) = — exp[-(v - E) jcr ] and p v (y) = = -, respectively. The 

cr Vtt ncr 2 + (v- E) 2 

magnitude of cr in the distribution function represents the strength of disorder. With p{6) in the 
hand, we can use Eq. (fT8l) to calculate the localization length numerically for cases of uncorrelated 
disorder. 

In our method, errors come from the value we choose as the number of discrete components 
of p{6) and where the cut-off is in the infinite range of the integration with respect to v. Basically, 
large E or small cr need large N, and large cr needs a large range of integration and a correction 
to remedy the effect of finite integration range. When is large, it spends a long time on the 

7 



(a) (b) 

FIG. 1. Probability density p{6) obtained numerically from Eq. (fTTT ) with E = (solid), 1 (dashed), 2 

1 1 

(dotted) when cr = 1 for (a) the Lorentzian distribution p v (v) = -, (b) the Gaussian distribution 

n 1 + (v - E) 1 

1 2 
p v (v) - — exp[-(v-£) ]. 

computation of p(0). In fact, this is always the main part of the computation time no matter how 
large N is. When the range of integration for v is large, more CPU time is consumed to compute 
the integration of v; increasing N does not improve the result when the range of integration remain 
unchanged; the result is sensitive to the range of integration of v for the Gaussian distribution, 
while for the Lorentzian distribution it is not. The reason should be that the Gaussian distribution 
varies more dramatically than the Lorentzian distribution. 



IV. RESULTS AND DISCUSSIONS 

As mentioned in the previous section, we obtained p(9) numerically from Eq. (fTTT) first. We 
performed the integral with respect to v by using Romberg's method, and solved the integral equa- 
tion by using a discrete matrix form, where p{6) was written in a vector form with N components. 
Figure \T\ shows the solutions of p{6) for the Lorentzian and Gaussian distributions respectively. 
When E = 0, p(9) is symmetric respected to = 0. As E becomes large, p{6) becomes a Dirac 
5-function. The position of the peak does not change monotonically with the increase of E. The 
salient difference of p{6) between the two distributions is, when E is small, the Lorentzian distri- 
bution has only one peak, while the Gaussian distribution has two. 

After we had the probability density p(8), the inverse localization length y was calculated from 
Eq. (fl~8"T) numerically. For the Lorentzian distribution, the analytical results were already available 



V(2 + Ef + cr 2 + V(2 - E) 2 + a 2 
y(E, cr) = arccosh . (23) 

Our numerical results for the Lorentzian distribution are shown in Figure|2l where we plot y versus 
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(a) (b) (c) 

FIG. 2. Comparison between the analytical and our numerical results for the Lorentzian distribution with cr 
- (a) 0.1, (b) 1, and (c) 10, respectively. 



FIG. 3. Inverse of the localization length y(E) obtained from Eq. (1201 and Eq. (1181 ) for the Gaussian (solid) 
and the Lorentzian (dashed) distributions. 

E with cr = 0.1, 1, 10 respectively with the comparison between our results and the analytical 
results. We found that they all coincide in excellence, and we can not see the difference from the 
figure. We used the Romberg integration method to perform a numerical integral with respect to v, 
and we made a cut-off to the infinity range of integration, so the real range is finite from E - 10 3 cr 
to E + 10 3 <x. Because the Lorentzian distribution does not decay rapidly to zero when cr was 
large, we need to consider the compensation of the contribution from the cut-out range. In fact, 
the integral of the cut-out range can be approximately integrated analytically, and the correction 
was 4 In |10 3 <xcos 0|/lO 3 7r. This correction would eliminate a constant difference between the 



numerical result and the analytical result. For weak disorder cr = 0.1 shown in Figure [2(a)} with 
increasing E, the number of mesh points N should go from 2000 to 9000 in order to obtain the right 
result, otherwise the result y will be much smaller than the analytical result. For strong disorder 



cr = 10 shown in Figure 2(c) , the number of mesh points N = 2000 with the correction is good 
enough to reach the desired accuracy, no matter how large E is. 

Numerical calculation of Eq. (1201) and Eq. (fl~8l) gave the same result in Figure [3] This is one 
evidence of the conclusion that z and tan 8 are the same quantity when the chain is long enough 
in the localization regime. From Figure [31 it can be seen that when in the band E = to 2, 
the difference of y between the Lorentzian and the Gaussian distributions is nearly a constant. 



(a) (b) 

FIG. 4. Inverse of the localization length y(cr) from Eq. (fT8l > for (a) the Lorentzian distribution and (b) the 
Gaussian distribution with E = 0, 1, 2, 3, 4 (from bottom to top). 



It seems that they have the same behavior in the band. When E becomes large, the y of the 
two distributions increase and go to the same value. The reason is that when E — » oo, both 
the Lorentzian and Gaussian distributions become a ^-function. For example, for the Lorentzian 
distribution, let y = v/E - 1, 

l im -7 — e?t— ~ dv = 6( y^ (24) 

£^c» ^ (v - £) 2 + cr 2 

then 

/(v)p v (v)dv = I f(Ey + E)6(y)dy 

oo *J — oo 

= /(£)■ (25) 

Thus if p{6) at large £ are the same for different p v (v), y{E) will be the same. We have checked 
that this is true for the two distributions. 

In Figure |4] we plot y versus cr. For the Lorentzian distribution, the results also coincide well 
with Eq. (|23T ). And we found that all the five lines can be fitted well by y{&) = a + b y/cr + ccr + 
do 2 . For all the five lines, d is always at least one order smaller than b and c. For the Gaussian 
distribution, this is also true when E is in the band. When E > 2 out of the band, the behavior 
of the Gaussian distribution is very different and it can not be fitted by the above formula. When 



cr becomes large, y(cr) becomes the same for different E. This is apparent in Figure [4(a)1 for the 
Lorentzian distribution. Similar to the derivation of large E, when cr grows large, the difference 
between p v (v) with different E vanishes. And this can be understood that when the width of the 
distribution becomes large, the position of the peak is irrelevant. 



In Il26ll the authors provided analytical results for the Gaussian distribution with weak and 
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E 


Lorentzian 


Gaussian 





0.49(r 


0.056cr 2 


1 


0.55(r 


0.08 lcr 2 


2 


0.72 


0.24cr 2/3 


3 


0.96 + 0.1 lo- 2 


0.96 - 0.063(r 2 


4 


1.32 + 0.045cr 2 


1.32 -0.022c- 2 



TABLE I. Fit for weak disorder < cr < 1 
strong disorder, which are summarized as follows, 



y(E, cr) = 



W 2 / 105.045- •• 
W 2 /24(4 - E 2 ) 
0.289- --iS 2 ) 1 ' 3 



E = 

< E < 2 . 
E = 2 



(26) 



In [|26[], the strength of disorder are represented by parameters W or 6, which can be expressed in 
terms of cr in our paper by W 2 = 6cr 2 and S 2 = cr 2 /2. 

For the purpose of comparison, we showed the fitting formulae of our numerical data at weak 
disorder < cr < 1 in Table [fl The ratio of the coefficients of E = 1 and E = derived from 
Eq. (|26l) is approximately 1 .46 and the ratio of our numerical results calculated from Table U is 
approximately 1.45. At the band edge E = 2, our result is y(2) ~ 0.24<x 2/3 , i.e. approximately 
0.30o" 2/3 . Therefore, our results confirmed the anomalous behavior of weak disorder at the band 



center E = |27 



29l] and the band edge E = 2 12811 . Furthermore, our numerical result gave an 



expression for strong disorder at E = 1 , y(cr) = -0.80 + 1 .00 In cr, which is in excellent agreement 

with the result of Il26n . where y(cr) = -0.797 h In cr when E is in the band. 

We noticed that for the Lorentzian distribution, the results can be obtained from the exact 
expression Eq. (1231) for small and large cr limits respectively. For strong disorder, y(cr) = In cr 
and our fitting formulae is y(cr) = 1 .00 In cr. It seems that the Lorentzian distribution should not 
have the anomalous behavior at the band center or band edge. However, it is apparent from Table 
H that at E = 2, the Lorentzian distribution also has an anomalous behavior similar to that of 
the Gaussian distribution. On the other hand, it should be mentioned that E = may should be 
seen as a boundary of two bands rather than the center of one band [11]. As shown in Figure [D 
the Lorentzian distribution has only one peak, while the Gaussian distribution has two when E is 
small. 
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V. CONCLUSIONS 



In this paper, we derived a parametrization method to deal with the transfer matrix of the one- 
dimensional Anderson model with diagonal uncorrelated disorder. With this method, we directly 
calculated the localization length under the thermodynamic limit in the localization regime. It 
avoids the difficulties faced by the traditional transfer matrix method; and without the sampling 
process, the accuracy can be improved easily. As we showed, the results of our method coincide 
very well with the known analytical results of the Lorentzian and the Gaussian distributions, in- 
cluding the anomalous behaviors at the band center and the band edge. It is quite efficient when 
the distribution of diagonal disorder is nonsingular, especially for moderate disorder. Further- 
more, we found that the Lorentzian distribution should give clues to the anomalies in the Gaussian 
distribution. Although it faces some difficulties for the cases like off-diagonal disorder or cor- 
related disorder, this method can be generalized to the coupled multichain system with diagonal 
uncorrelated disorder. 
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